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Abstract 

The frozen domain Effective Fragment Molecular Orbital method {PLoS ONE, 
accepted) is extended to allow for the treatment of a single fragment at the MP2 level 
of theory. The approach is applied to the conversion of chorismate to phrephenate 
by chorismate mutase, where the substrate is treated at the MP2 level of theory 
while the rest of the system is treated at the RHF level. MP2 geometry optimization 
is found to lower the barrier by up to 3.5 kcal/mol compared to RHF optimzations 
and ONIOM energy refinement and leads to smoother convergence with respect to 
basis set for the reaction profile. For double zeta basis sets the increase in CPU 
time relative to RHF is roughly a factor of two. 

1 Introduction 

Combined quantum mechanical/molecular mechanical and fragment-based quantum me- 
chanical methods^ are becoming increasingly popular in calculating properties of large 



2 



molecular structures. Recently, the Effective Fragment Molecular Orbital Method (EFMO)'^ 
was used to efficiently map the reaction path of the conversion of chorismate to prephanate 
in chorismate mutase. Via an additional ONIOM calculation step this reaction path was 
used to calculate the reaction barrier at a correlated (MP2) level of theory.l^ 

Such mappings require a stable and efficient algorithm for optimizing large structures 
in order to ensure convergence. The Frozen Domain (FD) and the Frozen Domain and 
Dimers (FDD) approximations are two multi-layer schemes, in which parts of the geometry 
are kept fixed during geometry optimization while the geometry of other parts are being 
optimized.'^ This reduces the CPU cost of the geometry optimization as well as it increases 
the optimization convergence since the search space is greatly reduced. 

In this work, we extend EFMO to work in a layer scheme in the FD or FDD ap- 
proximation, where the top layer consists of a single fragment (for instance a fragment 
containing a reaction complex) for which the one-body gradient contribution is calcu- 
lated at a correlated (MP2) level of theory, while the remaining one-body terms gradient 
terms are calculated using in a non-correlated manner (RHF). The presented method is 
discussed in detail in the next section. 

2 Theory 

The EFMO method partitions the wavefunction of the molecular system into fragments 
for which the wavefunction can be calculated separately or approximated by the classical 
Effective Fragment Potential. ^ ^ Geometry optimizations in this work are performed using 
the FDD approximations.!^ We note that the FD and FDD approximations and the use 
of EFMO in the FDD approximation is discussed in great detail elsewhere,'^ and only 
differences between the conventional implementation of FD/EFMO and our extension are 
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thus highlighted in the following. For a given molecular system, we define two layers F 
and A. The layer F is defined as all atoms having a frozen geometry and the layer A 
is defined as all atoms with a non-frozen geometry. Each layer is further divided into a 
number of molecular fragments. In this case, the EFMO energy is given as 



^EFMo ^ + + + (i) 

where E^ and E\ are the internal energies of region F and A, respectively, -Ep/^ is the 
interregion interaction between layers F and A, and E^^^ is the classical total polarization 
energy between fragments. In our EFM0-RHF:MP2 extension, we evaluate the internal 
energies of layers F and A at the RHF level. Furthermore, we specify a single fragment, 
S, (from the active region) to be treated at the MP2 level of theory. The total EFMO- 
RHF:MP2 energy is then given as 



77'EFMO-RHF:MP2 _ 7-.0,RHF , pO.RHF , pO.RHF , jt.POL , 77'0,MP2 /r,\ 
— £/p -h -h -C/p/A -^tot -^SeA ' l^J 



where E^^^"^ is the MP2 correlation energy of fragment S. 



The corresponding EFMO energy gradients of each region in the FD approximation: 

(3) 



dE™'' _ dEi ggj/p a£;poL 

dxA dxA dx\ dxx 
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Table 1: Electronic energy barrier for chorismate mutase and the corresponding reaction 
coordinate for the transition state. EFM0-RHF/MP2 results are from the presented 
work, while ONIOM results are obtained from Steinmann et al^ 
Method i?(TS) Energy barrier Reaction energy 



EFMO-RHF/6-31G(d):MP2/6-31G(d) 

EFMO-RHF /6-31G (d) :MP2 /cc-pVDZ 

EFMO-RHF /6-31G (d) :MP2 /cc-pVTZ 

ONIOM(MP2/6-31G(d):RHF/6-31G(d)ys 

ONIOM(MP2/cc-pVDZ:RHF/6-31G(d))ISl 

ONIOM(MP2/cc-pVTZ:RHF/6-31G(d))E' 

ONIOM(MP2/cc-pVQZ:RHF/6-31G(d))S 



-0.17 A 20.95 kcal/mol -4.79 kcal/mol 

-0.43 A 19.21 kcal/mol -6.83 kcal/mol 

-0.43 A 18.34 kcal/mol -6.17 kcal/mol 

0.13 A 22.24 kcal/mol -3.20 kcal/mol 

-0.36 A 19.75 kcal/mol -5.48 kcal/mol 

0.13 A 21.79 kcal/mol -1.17 kcal/mol 

0.13 A 21.68 kcal/mol -0.82 kcal/mol 



This gives the following EFM0-RHF:MP2 energy gradients: 



^^EFMO-RHF:MP2 ^^EFMO ^^MP2 

s = — + ^r^ 

oxa oxa ox a 

a^EFMO-RHF:MP2 n^EFMO 

= (6) 



dxv dxi 



Where ^^^^ contains the MP2 part of the gradient for fragment SgA (and not the 
Hartree-Fock part). 



3 Results 

Reaction Barrier. Electronic energy barriers and reaction coordinates for the transition 
state are given in Table [3] and Fig. [3] We find the electronic energy barrier at the EFMO- 
RHF/6-31G(d):MP2/6-31G(d) level of theory to be 20.95 kcal/mol. Increasing the size 
of the basis set on the MP2 fragment decreases the barrier to 19.21 kcal/mol with the 
cc-pVDZ basis set and 18.34 kcal/mol with the cc-pVTZ basis set. 

In comparison, the corresponding MP2:RHF ONIOM calculations by Steinmann et 
al. resulted in barriers of 22.24, 19.75, and 21.79 kcal/mol, respectively. In contrast to 
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Figure 1: F denotes the frozen geometry region (green); A denotes the active geometry 
region (red); S & A denotes fragment S, for which the MP2 energy and gradients are 
evaluated (yellow). 



6 



the ONIOM approach, we find that for increasing basis set sizes, the electronic energy 
barrier is systematically reduced. The experimental enthalpy barrier has been measured 
to be 12.7 kcal/mol.EEl 

Transition State Structure. We define the reaction coordinate similarly to Claeyssens 
et al^as the difference in bond length between the breaking 02-Cl bond and the forming 
C4-C3 bond in chorismate, i.e. R = R21 — -R43 (see Figjs]). The reaction coordinate of the 
transition state was found to be -0.17 A using the 6-31G(d) basis set on the MP2 fragment 
and -0.43 A for both the cc-pVTZ and cc-pVDZ basis set reaction paths. This conver- 
gence with respect to basis set is in good, quantitative agreement with the coordinates 
obtained obtained by Szefczyk et al^ 

In comparison, the corresponding MP2:RHF ONIOM calculations by Steinmann et 
a/.'^ resulted in transition state reaction coordinates of 0.13, -0.36, and 0.13 A with the 
cc-pVDZ, cc-pVTZ and cc-pVQZ basis sets used in the MP2 calculation, respectively 

Reaction Energy. The energy difference between the product and reactant state 
is found to be -3.2 kcal/mol using the 6-31G(d) basis set on chorismate. Increasing the 
basis set to cc-pVDZ and cc-pVTZ on chorismate decreased the reaction energy to -6.83 
kcal/mol and -6.17 kcal/mol, respectively. The ONIOM approach by Steinmann et al. 
found the reaction energy to be between -5.48 kcal/mol to -0.82 kcal/mol. However, in 
the ONIOM approach increasing the basis set from cc-pVTZ on chorismate increased the 
reaction energy from -5.48 kcal/mol to -1.17 kcal/mol. We find that all three basis sets are 
in close agreement, and only a 0.7 kcal/mol difference between the cc-pVDZ and cc-pVTZ 
reaction paths. 

Timings Running on 80 cores distributed on 10 compute ndes and using the de- 
fault compute node load balancing scheme, the average time for a geometry optimization 
step was 760 s at the EFMO-RHF/6-31G(d) level of theory. ^ For the EFMO-RHF/6- 
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Table 2: Timings for the average geometry optimization step for Chorismate mutase using 
using different metliods. EFMO-RHF/6-31G(d) timings are obtained from Steinmann et 
al^ The timing marked (distributed) denotes that in this calculation, the MP2 part was 



distributed across all nodes (see text). 
Method Average step time 

EFMO-RHF/6-31G(d):MP2/6-31G(d) 1527i 

EFMO-RHF/6-31G(d):MP2/cc-pVDZ 1967 s 

EFMO-RHF/6-31G(d):MP2/cc-pVTZ 18845 s 

EFMO-RHF/6-31G(d):MP2/cc-pVTZ (distributed) 10911 s 

EFMO-RHF/6-31G(d)E! 760 s 




1 2 3 



Figure 2: Claisen rearrangement of chorismate to prephenate. The atoms describing the 
reaction coordinate are marked with numbers one trough four. 



31G(d):MP2/6-31G(d) calculation, this time increased to 1526 s per step. Increasing the 
basis set on the MP2 part of the system to cc-pVDZ and cc-pVTZ increased the time to 
1967 s and 18845 s, respectively (see Table |3]). The large increase in calculation time from 
cc-pVDZ to cc-pVTZ was found to be due to sub-optimal load balancing during the MP2 
part of the calculation. Subsequently, one optimization using the cc-pVTZ was carried 
out, in which the calculation of the MP 2 fragment energy and gradient was distributed 
across all 10 nodes. This reduced the average gradient step time from 18845 s to 10911 s. 
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EFMO-RHF/6-31G{d):MP2/? electronic energy barrier 

25 1 \ 1 \ \ 1 \ 1 




_i O' ' I ' ' I ' I I 

-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 

Reaction coordinate 



Figure 3: Electronic energy versus reaction coordinate for the convesion of chorismate 
to prephanate in chorismate mutase. The three reation paths are calculated using the 
FDD/EFM0-RHF:MP2 approach with three different basis sets on the reaction complex 
in the MP2 region. The 6-31G(d) basis set was used on the RHF regions in all three cases. 
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4 Methods 

All calculations were carried out in a development version of GAMESS'^ containing ver- 
sion 4.3 of the EFMO module. 

Starting structures for Chorismate Mutase were obtained from Steinmann et alJ^ who 
prepared the structures following Claeyssens et alJ^ The structures were prepared for 
fragment based calculations in Fraglt.'^^' The adiabatic mapping was carried out using 
the presented EFM0-RHF:MP2 gradient with 6-31G(d) basis set on all atoms. Two 
additional runs were also carried out, in these cases with the cc-pVDZ or cc-pVTZ on 
chorismate and 6-31G(d) on remaining atoms. 

Points on the EFMO-RHF/6-31G(d):MP2//cc-pVTZ reaction path were started from 
the converged structures with identical reaction coordinates in the EFMO-RHF/6-31G(d):MP2/ /cc- 
pVDZ reaction path. All residues with an atom within a distance of 2.0 A from any atom 
in chorismate were assigned to the A (active) region. 

The RESDIM keyword was set to 1.5 and the optimization convergence criterion was 
set to 5.0 ■ 10~^ Hartree/Bohr. Each step of the reaction path was obtained by imposing 
harmonic constraints on R12 and -R13 with a force constant of 500 kcal/A. The FDD 
approximation was enabled by setting M0DFD=3 and was enabled for all calculations. 
Our code is enabled by specifying MPLEVL(1)=0,2 in in the $FMO keyword group and 
using a new keyword ISFRG=A^ in the $FMO keyword group which specifies that the 
one-body properties of fragment N are to be evaluated at the MP2 level, while one-body 
properties of other fragments in the A layer will be calculated at the RHF level. 

Timings for the optimization procedure were carried out on 80 Intel Xeon X5550 CPU 
cores distributed across 10 nodes. 
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5 Conclusion 

We have implemented an scheme for optimizing a reaction complex using a correlated 
method in the FD /EFMO approximation. ^ Our method is computationally efficient when 
a moderately sized basis sets is used on the correlated fragment or if the calculation is 
not distributed across different compute nodes. 

While our EFM0-RHF:MP2 approach does not achieve chemical accuracy in predict- 
ing enthalpy barrier of the conversion of chorismate to prephanate in chorismate mutase, 
we have demonstrated that our method serves as a rigorous and viable alternative to the 
widely used ONIOM approach. 
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